%% Fig.1 STEP SCHEMATIC
%(a) initial bathymetry
figure('color','white')
set(0,'DefaultAxesFontName','Calibri','DefaultAxesFontSize',7,...
    'DefaultLineLineWidth',1);
set(gcf,'Units','centimeters','Position',[0 0 6 3]);
colormap(colormap_cpt('ETOPO1')); 
k = [-10 10];
contourf(X(1:300,1:300)/1000,Y(1:300,1:300)/1000,ZI{1}(1:300,1:300),100,'linecolor','none');caxis(k);daspect([1 1 1]); hold on; axis on;axis('off');
print('-dtiff','-r300','initial bathymetry')

%(b) equilibrium bathymetry
figure('color','white')
set(gcf,'Units','centimeters','Position',[0 0 6 6]);
colormap(colormap_cpt('ETOPO1')); 
k = [-10 10];colorbar('location','southoutside','Fontname','Times','Fontsize',7);
contourf(X(1:300,1:300)/1000,Y(1:300,1:300)/1000,Z1{1}(1:300,1:300),100,'linecolor','none');caxis(k);daspect([1 1 1]); hold on; axis on;axis('off');colorbar('location','southoutside','Fontname','Times','Fontsize',7);
print('-dpng','-r300','initial bathymetry2')

%(c) SLR value
figure('color','white')
set(gcf,'Units','centimeters','Position',[0 0 4 3]);
plot((0:10:100),SLR(:,2),'linewidth',1,'color',[1 1 1]*0.5);
hold on;
plot((0:10:100),SLR(:,3),'-.','linewidth',1,'color',[1 1 1]*0.5);
plot((0:10:100),SLR(:,4),'linewidth',1,'color',[1 1 1]*0);
plot((0:10:100),SLR(:,5),'-.','linewidth',1,'color',[1 1 1]*0);
plot((0:10:100),zeros(11,1),'-.','linewidth',1,'color',[1 1 1]*0.7); grid on;
ylim([-0.2 1]);
leg = legend('Scenario 1','Scenario 2','Scenario 3','Scenario 4','Base','location','northwest','Fontsize',5,'box','off');
leg.ItemTokenSize=[10,3];
xlabel('Time (year)','fontsize',7); ylabel('SLR (m)','fontsize',7)
set(gca,'Fontsize',7,'Fontname','Times')
saveas(gcf,[figsavedir 'SLR value' '.fig']);
print('-dtiff','-r600','SLR value')

%(d) flood curve
figure('color','white')
set(gcf,'Units','centimeters','Position',[0 0 4 3]);
plot((floodcurve(:,1)-floodcurve(2,1))/1440,floodcurve(:,2)*6,'k');grid on;
xlabel('Time (day)','fontsize',7); ylabel('Discharge (m^3/s)','fontsize',7);
set(gca,'Fontsize',7,'Fontname','Times')
ylim([0 3000])
print('-dpng','-r600','flood curve')

%% Fig.2 model validation
%calculation
for i = 1:20
    Prism0(i) = sum(Q1{i}(1:149).*(Q1{i}(1:149)<0))*-5*60;
    CSA0(i) = sum(Z1{1}(1:300,i*10).*(Z1{1}(1:300,i*10)<0))*-50;
end

figure('color','white')
set(gcf,'Units','centimeters','Position',[0 0 15 5]);
subplot(1,2,2)
plot(0.1:0.1:20,smooth(mean(ZI{1}(1:300,1:200).*active1(1:300,1:200),'omitnan')),'linewidth',1.5); hold on;
plot(0.1:0.1:19.7,smooth(mean(ZI{22}(1:300,1:197).*active1(1:300,1:197),'omitnan')),'linewidth',1.5); hold on;
plot(0.1:0.1:19.7,smooth(mean(ZI{45}(1:300,1:197).*active1(1:300,1:197),'omitnan')),'linewidth',1.5);
plot(0.1:0.1:19.7,smooth(mean(ZI{67}(1:300,1:197).*active1(1:300,1:197),'omitnan')),'linewidth',1.5);
plot(0.1:0.1:19.7,smooth(mean(ZI{89}(1:300,1:197).*active1(1:300,1:197),'omitnan')),'linewidth',1.5);
ylabel('Bed elevation (m)'); xlabel('Distance to upstream boundary (km)')
set(gca,'Fontsize',7,'Fontname','Times'); grid on;
legend('250 year','500 year','750 year','1000 year','location','Southwest','Fontsize',7)
subplot(1,2,1)
scatter(CSA0(:),Prism0(:),'filled'); hold on;
set(gca,'xscale','log');set(gca,'yscale','log');
plot(10:1:1000,(10:1:1000)/(1.576*10^-4)^(1/0.95),'--','color',[.2 .2 .2],'linewidth',1.5)
plot(10:1:1000,(10:1:1000)/(3.039*10^-5)^(1/1.05),'-.','color',[.2 .2 .2],'linewidth',1.5)
plot(10:1:1000,(10:1:1000)/(9.311*10^-4)^(1/0.84),'--','color',[.7 .7 .7],'linewidth',1.5)
plot(10:1:1000,(10:1:1000)/(2.833*10^-4)^(1/0.91),'-.','color',[.7 .7 .7],'linewidth',1.5)
plot(10:1:1000,(10:1:1000)/(6.56*10^-5)^(1/1),':','color',[.2 .2 .2],'linewidth',1.5)
set(gca,'Fontsize',7,'Fontname','Times');
h = legend('Model result','Jarrett (1976)(all inlets)','Jarrett (1976)(Atlantic coast)','Jarrett (1976)(Gulf coast)','Jarrett (1976)(Pacific coast)','O Brien (1969)','location','Southeast','Fontsize',5);
h.ItemTokenSize=[10,3];
xlabel('Cross-section area (m^2)');ylabel('Tidal prism (m^3)');
ylim([100 10^8])
annotation('textbox',...
    [0.6 0.83 0.042 0.086],...
    'String',{'(b)'},...
    'FontName','times',...
    'EdgeColor','none');
annotation('textbox',...
    [0.15 0.83 0.042 0.086],...
    'String',{'(a)'},...
    'FontName','times',...
    'EdgeColor','none');
 print('fig2.eps','-depsc2');

%% Fig.3 flood_water_level difference
figure('color','white')
x= 0.1; y = 0.9; dx = 0.167; dy = 0.215;
annotation('textbox',[x+0*dx y 0.06 0.05],'String',{'(a)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+1*dx y 0.06 0.05],'String',{'(b)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+2*dx y 0.06 0.05],'String',{'(c)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+3*dx y 0.06 0.05],'String',{'(d)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+4*dx y 0.06 0.05],'String',{'(e)'},'EdgeColor','none','Fontname','Times','Fontsize',8);

annotation('textbox',[x y-1*dy 0.06 0.05],'String',{'(f)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x y-2*dy 0.06 0.05],'String',{'(k)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x y-3*dy 0.06 0.05],'String',{'(p)'},'EdgeColor','none','Fontname','Times','Fontsize',8);

annotation('textbox',[x+dx y-1*dy 0.06 0.05],'String',{'(g)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+dx y-2*dy 0.06 0.05],'String',{'(l)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+dx y-3*dy 0.06 0.05],'String',{'(q)'},'EdgeColor','none','Fontname','Times','Fontsize',8);

annotation('textbox',[x+2*dx y-1*dy 0.06 0.05],'String',{'(h)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+2*dx y-2*dy 0.06 0.05],'String',{'(m)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+2*dx y-3*dy 0.06 0.05],'String',{'(r)'},'EdgeColor','none','Fontname','Times','Fontsize',8);

annotation('textbox',[x+3*dx y-1*dy 0.06 0.05],'String',{'(i)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+3*dx y-2*dy 0.06 0.05],'String',{'(n)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+3*dx y-3*dy 0.06 0.05],'String',{'(s)'},'EdgeColor','none','Fontname','Times','Fontsize',8);

annotation('textbox',[x+4*dx y-1*dy 0.06 0.05],'String',{'(j)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+4*dx y-2*dy 0.06 0.05],'String',{'(o)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x+4*dx y-3*dy 0.06 0.05],'String',{'(t)'},'EdgeColor','none','Fontname','Times','Fontsize',8);

annotation('textbox',[0 y-0*dy-0.18 0.18 0.18],'String',{'Scenario 1'},'EdgeColor','none','Fontname','Times');
annotation('textbox',[0 y-1*dy-0.18 0.18 0.18],'String',{'Scenario 2'},'EdgeColor','none','Fontname','Times');
annotation('textbox',[0 y-2*dy-0.18 0.18 0.18],'String',{'Scenario 3'},'EdgeColor','none','Fontname','Times');
annotation('textbox',[0 y-3*dy-0.18 0.18 0.18],'String',{'Scenario 4'},'EdgeColor','none','Fontname','Times');

set(gcf,'Units','centimeters','Position',[0 0 18.0 7.5]);
colormap(colormap_cpt('BlWhRe'));
subplot(4,5,5)
k = [-0.1 0.1];%title('Scenario1')
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif1(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
title('100 year');
set(gca,'xtick',[])
set(gca,'ytick',[])
subplot(4,5,10)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif2(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'xtick',[])
set(gca,'ytick',[])
subplot(4,5,15)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif3(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'xtick',[])
set(gca,'ytick',[])
a = colorbar('position',[.91 .1 .015 .85]);
a.Label.String = 'Flood water level difference (m)';
a.Label.Position = [3.2 0];
subplot(4,5,20)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif4(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'ytick',[])
subplot(4,5,1)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif2012(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
title('20 year');%ylabel('Scenario 1');
set(gca,'xtick',[])
subplot(4,5,6)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif2012(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'xtick',[])
%ylabel('Scenario 2');
subplot(4,5,11)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif2034(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'xtick',[])
%ylabel('Scenario 3');
subplot(4,5,16)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif2034(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
%ylabel('Scenario 4');
subplot(4,5,2)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif4012(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
title('40 year');
set(gca,'xtick',[])
set(gca,'ytick',[])
subplot(4,5,7)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif4012(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'xtick',[])
set(gca,'ytick',[])
subplot(4,5,12)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif4034(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'xtick',[])
set(gca,'ytick',[])
subplot(4,5,17)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif4034(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'ytick',[])
subplot(4,5,3)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif601(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
title('60 year');
set(gca,'xtick',[])
set(gca,'ytick',[])
subplot(4,5,8)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif602(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'xtick',[])
set(gca,'ytick',[])
subplot(4,5,13)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif603(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'xtick',[])
set(gca,'ytick',[])
subplot(4,5,18)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif604(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'ytick',[])
xlabel('Distance (km)');
subplot(4,5,4)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif801(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
title('80 year');
set(gca,'xtick',[])
set(gca,'ytick',[])
subplot(4,5,9)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif802(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'xtick',[])
set(gca,'ytick',[])
subplot(4,5,14)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif803(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'xtick',[])
set(gca,'ytick',[])
subplot(4,5,19)
contourf(X(75:225,1:220)/1000,Y(75:225,1:220)/1000,wl_dif804(75:225,1:220),300,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontsize',7,'Fontname','Times');
set(gca,'ytick',[])
subaxis(4, 6, 1, 'sh', 0.03, 'sv', 0.01, 'padding', 0, 'margin', 0);
print fig3.eps -depsc2 -r600
print('-dpng','-r300','fig3')

%% Fig 4 flood inundation

figure
set(gcf,'Units','centimeters','position',[0 0 7 5])
plot((1:1:151)/4*745/24/60,inundate1,'linewidth',1.1); hold on;
plot((1:1:151)/4*745/24/60,inundate2,'linewidth',1.1);
plot((1:1:151)/4*745/24/60,inundate3,'linewidth',1.1);
plot((1:1:151)/4*745/24/60,inundate4,'linewidth',1.1);
plot((1:1:151)/4*745/24/60,inundate5,'linewidth',1.1); ylim([0 4*10^7]);
ylabel('Inundation Area (m^2)'); xlabel('Time (day)')
legend('Base','Scenario 1','Scenario 2','Scenario 3','Scenario 4')
set(gca,'Fontsize',7,'Fontname','Times');
print('fig4.eps','-depsc2');

%% Fig 5 water level decompse
figure('color','white')
set(gcf,'Units','centimeters','Position',[0 0 14.0 16.0]);

x= 0.13; y = 0.905; dx = 0.2; dy = 0.22;
annotation('textbox',[x y 0.06 0.05],'String',{'(a)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x y-dy 0.06 0.05],'String',{'(b)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x y-dy*2 0.06 0.05],'String',{'(c)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x y-dy*3 0.06 0.05],'String',{'(d)'},'EdgeColor','none','Fontname','Times','Fontsize',8);

k = [-0.3 0.3];
subplot(4,1,4);hold on;
title('scenario 4')
bar(1:20,[((2500./C5./Wid1./S1.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid5./S1.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid1./S5.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte5)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1)]);
 plot(((2500./C5./Wid5./S5.^0.5).^(2/3)+Bte5)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1),'-o','markersize',4,'linewidth',1.5);ylim(k);grid on;
 plot(wld5-wld1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'ytick',-0.3:0.1:0.3); set(gca,'Fontname','Times','Fontsize',7);
xlabel('Distance to upstream boundary (km)'); ylabel('water level difference (m)');
subplot(4,1,3); hold on;
title('scenario 3')
bar(1:20,[((2500./C4./Wid1./S1.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid4./S1.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid1./S4.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte4)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1)]);
 plot(((2500./C4./Wid4./S4.^0.5).^(2/3)+Bte4)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1),'-o','markersize',4,'linewidth',1.5);ylim(k);grid on;
 plot(wld4-wld1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'ytick',-0.3:0.1:0.3); set(gca,'Fontname','Times','Fontsize',7);
 ylabel('water level difference (m)');
subplot(4,1,2); hold on;
title('scenario 2')
bar(1:20,[((2500./C3./Wid1./S1.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid3./S1.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid1./S3.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte3)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1)]);
 plot(((2500./C3./Wid3./S3.^0.5).^(2/3)+Bte3)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1),'-o','markersize',4,'linewidth',1.5);ylim(k);grid on;
 plot(wld3-wld1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'ytick',-0.3:0.1:0.3); set(gca,'Fontname','Times','Fontsize',7);
  ylabel('water level difference (m)');
subplot(4,1,1); hold on; 
title('scenario 1')
bar(1:20,[((2500./C2./Wid1./S1.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid2./S1.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid1./S2.^0.5).^(2/3)+Bte1)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1) 
                 ((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte2)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1)]);
 plot(((2500./C2./Wid2./S2.^0.5).^(2/3)+Bte2)-((2500./C1./Wid1./S1.^0.5).^(2/3)+Bte1),'-o','markersize',4,'linewidth',1.5);ylim(k);grid on;
  ylabel('water level difference (m)');
 plot(wld2-wld1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'ytick',-0.3:0.1:0.3);
 h = legend('Chezy Roughness','Channel Width','Surface Slope','Bed Elevation','Total difference','Model result','location',[0.26,0.87,0.55,0.05],'NumColumns',3,'Box','off');
 h.ItemTokenSize = [14,7];
 set(gca,'Fontname','Times','Fontsize',7);
 print('fig5.eps','-depsc2');

 print('-dsvg','fig5_water_level_decompose')
 
 
 figure('color','white')
set(gcf,'Units','centimeters','Position',[0 0 14.0 16.0]);
k = [-0.3 0.3];
subplot(4,1,4);hold on;
title('scenario 4')
bar(1:20,[((2500.*M5./Wid1./S1.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid5./S1.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid1./S5.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte5)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1)]);
 plot(((2500.*M5./Wid5./S5.^0.5).^(3/5)+Bte5)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1),'-o','markersize',4,'linewidth',1.5);ylim(k);grid on;
 plot(wld5-wld1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'ytick',-0.3:0.1:0.3); set(gca,'Fontname','Times','Fontsize',7);
xlabel('Distance to upstream boundary (km)'); ylabel('water level difference (m)');
subplot(4,1,3); hold on;
title('scenario 3')
bar(1:20,[((2500.*M4./Wid1./S1.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid4./S1.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid1./S4.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte4)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1)]);
 plot(((2500.*M4./Wid4./S4.^0.5).^(3/5)+Bte4)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1),'-o','markersize',4,'linewidth',1.5);ylim(k);grid on;
 plot(wld4-wld1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'ytick',-0.3:0.1:0.3); set(gca,'Fontname','Times','Fontsize',7);
 ylabel('water level difference (m)');
subplot(4,1,2); hold on;
title('scenario 2')
bar(1:20,[((2500.*M3./Wid1./S1.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid3./S1.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid1./S3.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte3)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1)]);
 plot(((2500.*M3./Wid3./S3.^0.5).^(3/5)+Bte3)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1),'-o','markersize',4,'linewidth',1.5);ylim(k);grid on;
 plot(wld3-wld1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'ytick',-0.3:0.1:0.3); set(gca,'Fontname','Times','Fontsize',7);
  ylabel('water level difference (m)');
subplot(4,1,1); hold on; 
title('scenario 1')
bar(1:20,[((2500.*M2./Wid1./S1.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid2./S1.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid1./S2.^0.5).^(3/5)+Bte1)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1) 
                 ((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte2)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1)]);
 plot(((2500.*M2./Wid2./S2.^0.5).^(3/5)+Bte2)-((2500.*M1./Wid1./S1.^0.5).^(3/5)+Bte1),'-o','markersize',4,'linewidth',1.5);ylim(k);grid on;
  ylabel('water level difference (m)');
 plot(wld2-wld1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'ytick',-0.3:0.1:0.3);
 h = legend("Manning's Roughness",'Channel Width','Surface Slope','Bed Elevation','Total difference','Model result','location',[0.26,0.87,0.55,0.05],'NumColumns',3,'Box','off');
 h.ItemTokenSize = [14,7];
 set(gca,'Fontname','Times','Fontsize',7);
 print('-dpng','-r300','fig5_water_level_decompose')
 
 %% Fig6 SLR influence on morphology
figure('color','white')
set(gcf,'Units','centimeters','Position',[0 0 15.0 10.0]);
%colormap(colormap_cpt('ETOPO1'));
k = [-10 10];
subplot(5,2,1);
contourf(X(100:200,2:200)/1000,Y(100:200,2:200)/1000,Z1{283}(100:200,2:200),'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontname','Times','Fontsize',7);
annotation('textbox',[0.08 0.925 0.06 0.05],'String',{'(a)'},'EdgeColor','none','Fontname','Times');
colormap(gca,colormap_cpt('ETOPO1')); ylabel('Base')
subplot(5,2,3);
contourf(X(100:200,2:200)/1000,Y(100:200,2:200)/1000,Z2{283}(100:200,2:200),'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontname','Times','Fontsize',7);
annotation('textbox',[0.08 0.765 0.06 0.05],'String',{'(b)'},'EdgeColor','none','Fontname','Times');
colormap(gca,colormap_cpt('ETOPO1'));ylabel('Scenario1')
subplot(5,2,5);
contourf(X(100:200,2:200)/1000,Y(100:200,2:200)/1000,Z3{283}(100:200,2:200),'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontname','Times','Fontsize',7);
annotation('textbox',[0.08 0.595 0.06 0.05],'String',{'(c)'},'EdgeColor','none','Fontname','Times');
colormap(gca,colormap_cpt('ETOPO1'));ylabel('Scenario2')
subplot(5,2,7);
contourf(X(100:200,2:200)/1000,Y(100:200,2:200)/1000,Z4{283}(100:200,2:200),'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontname','Times','Fontsize',7);
annotation('textbox',[0.08 0.425 0.06 0.05],'String',{'(d)'},'EdgeColor','none','Fontname','Times');
colormap(gca,colormap_cpt('ETOPO1'));ylabel('Scenario3')
subplot(5,2,9);
contourf(X(100:200,2:200)/1000,Y(100:200,2:200)/1000,Z5{283}(100:200,2:200),'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'Fontname','Times','Fontsize',7);
annotation('textbox',[0.08 0.255 0.06 0.05],'String',{'(e)'},'EdgeColor','none','Fontname','Times');
colormap(gca,colormap_cpt('ETOPO1'));xlabel('Distance (km)'); ylabel('Scenario4')
colorbar('position',[0.48,0.12,0.012,0.8])
k = [-1 1];%colormap(colormap_cpt('BlWhRe'));
subplot('position',[.585 .64 .3 .115 ]);
contourf(bte2'-bte1',100,'linecolor','none');caxis(k);title('��Z')
annotation('textbox',[0.54 0.765 0.06 0.05],'String',{'(f)'},'EdgeColor','none','Fontname','Times');
set(gca,'ytick',0:283/4:283,'yticklabel',0:25:100,'xtick',0:50:200,'xticklabel',0:5:20);colormap(gca,colormap_cpt('BlWhRe'));set(gca,'Fontname','Times','Fontsize',7);
subplot('position',[.585 .467 .3 .115 ]);
contourf(bte3'-bte1',100,'linecolor','none');caxis(k);
annotation('textbox',[0.54 0.595 0.06 0.05],'String',{'(g)'},'EdgeColor','none','Fontname','Times');
set(gca,'ytick',0:283/4:283,'yticklabel',0:25:100,'xtick',0:50:200,'xticklabel',0:5:20);colormap(gca,colormap_cpt('BlWhRe'));set(gca,'Fontname','Times','Fontsize',7);
subplot('position',[.585 .293 .3 .115 ]);
contourf(bte4'-bte1',100,'linecolor','none');caxis(k);
annotation('textbox',[0.54 0.425 0.06 0.05],'String',{'(h)'},'EdgeColor','none','Fontname','Times');
set(gca,'ytick',0:283/4:283,'yticklabel',0:25:100,'xtick',0:50:200,'xticklabel',0:5:20);colormap(gca,colormap_cpt('BlWhRe'));set(gca,'Fontname','Times','Fontsize',7);
subplot('position',[.585 .12 .3 .115 ]);
contourf(bte5'-bte1',100,'linecolor','none');caxis(k);xlabel('Distance (km)')
annotation('textbox',[0.54 0.255 0.06 0.05],'String',{'(i)'},'EdgeColor','none','Fontname','Times');
set(gca,'ytick',0:283/4:283,'yticklabel',0:25:100,'xtick',0:50:200,'xticklabel',0:5:20);colormap(gca,colormap_cpt('BlWhRe'));set(gca,'Fontname','Times','Fontsize',7);
colorbar('position',[0.92,0.12,0.012,0.61])
print fig6.jpg -djpeg -r600

figure('color','white')
set(gcf,'Units','centimeters','Position',[0 0 15.0 10.0]);
plot(Z1{:}(100:200,151))




%% Fig7 Tidal components
figure('color','white')
set(gcf,'Units','centimeters','Position',[0 0 20.0 20.0]);
subplot(3,1,2)
h1 = plot((0.1:0.1:19.9),mean(P1{2,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'linewidth',1.5,'color','#0072BD');
hold on;
h2 = plot((0.1:0.1:19.9),mean(P1{2,2}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'--','linewidth',1.5,'color','#0072BD');
h3 = plot((0.1:0.1:19.9),mean(P2{2,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'linewidth',1.5,'color','#D95319');
h4 = plot((0.1:0.1:19.9),mean(P2{2,2}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'--','linewidth',1.5,'color','#D95319');
h5 = plot((0.1:0.1:19.9),mean(P3{2,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'linewidth',1.5,'color','#EDB120');
h6 = plot((0.1:0.1:19.9),mean(P3{2,2}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'--','linewidth',1.5,'color','#EDB120');
h7 = plot((0.1:0.1:19.9),mean(P4{2,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'linewidth',1.5,'color','#7E2F8E');
h8 = plot((0.1:0.1:19.9),mean(P4{2,2}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'--','linewidth',1.5,'color','#7E2F8E');
h9 = plot((0.1:0.1:19.9),mean(P5{2,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'linewidth',1.5,'color','#77AC30');
h10 = plot((0.1:0.1:19.9),mean(P5{2,2}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'--','linewidth',1.5,'color','#77AC30');
hold on;
ylabel('Tidal amplitude (m)');set(gca,'Fontname','Times','Fontsize',10); grid on;
annotation('textbox',[0.12 0.91 0.06 0.05],'String',{'(a)'},'EdgeColor','none','Fontname','Times');
legend([h1 h3 h5 h7 h9 h2 h4 h6 h8 h10],{'M2 (Base)','M2(Sce 1)','M2 (Sce 2)','M2 (Sce 3)','M2 (Sce 4)','M4 (Base)','M4 (Sce 1)','M4 (Sce 2)','M4 (Sce 3)','M4 (Sce 4)'},'location','west','NumColumns',1);

subplot(3,1,3)
h1 = plot((0.1:0.1:19.9),mean(P1{7,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'linewidth',1.5,'color','#0072BD');
hold on;
h2 = plot((0.1:0.1:19.9),mean(P1{9,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'--','linewidth',1.5,'color','#0072BD');
h3 = plot((0.1:0.1:19.9),mean(P2{7,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'linewidth',1.5,'color','#D95319');
h4 = plot((0.1:0.1:19.9),mean(P2{9,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'--','linewidth',1.5,'color','#D95319');
h5 = plot((0.1:0.1:19.9),mean(P3{7,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'linewidth',1.5,'color','#EDB120');
h6 = plot((0.1:0.1:19.9),mean(P3{9,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'--','linewidth',1.5,'color','#EDB120');
h7 = plot((0.1:0.1:19.9),mean(P4{7,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'linewidth',1.5,'color','#7E2F8E');
h8 = plot((0.1:0.1:19.9),mean(P4{9,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'--','linewidth',1.5,'color','#7E2F8E');
h9 = plot((0.1:0.1:19.9),mean(P5{7,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'linewidth',1.5,'color','#77AC30');
h10 = plot((0.1:0.1:19.9),mean(P5{9,1}(1:300,1:199).*active1(1:300,1:199),'omitnan'),'--','linewidth',1.5,'color','#77AC30');
xlabel('Distance to upstream boundary (km)');
ylabel({'Velocity amplitude (m/s)'}); ylim([0 0.5]);set(gca,'Fontname','Times','Fontsize',10);grid on;
annotation('textbox',[0.12 0.61 0.06 0.05],'String',{'(b)'},'EdgeColor','none','Fontname','Times');


subplot(3,1,1)
plot((0.1:0.1:19.9),smooth(mean(P1{1,1}(1:300,1:199).*active2(1:300,1:199),'omitnan')),'-','linewidth',1.5,'color','#0072BD');hold on;
plot((0.1:0.1:19.9),smooth(mean(P2{1,1}(1:300,1:199).*active2(1:300,1:199),'omitnan')),'-','linewidth',1.5,'color','#D95319');
plot((0.1:0.1:19.9),smooth(mean(P3{1,1}(1:300,1:199).*active2(1:300,1:199),'omitnan')),'-','linewidth',1.5,'color','#EDB120');
plot((0.1:0.1:19.9),smooth(mean(P4{1,1}(1:300,1:199).*active2(1:300,1:199),'omitnan')),'-','linewidth',1.5,'color','#7E2F8E');
plot((0.1:0.1:19.9),smooth(mean(P5{1,1}(1:300,1:199).*active2(1:300,1:199),'omitnan')),'-','linewidth',1.5,'color','#77AC30');
ylim([0 0.6]);
ylabel('Mean flow velocity(m/s)')
set(gca,'Fontname','Times','Fontsize',10);
legend('Base','Scenario 1','Scenario 2','Scenario 3','Scenario 4','location','Southwest');grid on;
annotation('textbox',[0.12 0.31 0.06 0.05],'String',{'(c)'},'EdgeColor','none','Fontname','Times');
print fig7.eps -depsc2 -r600



figure
plot((0.1:0.1:19.9),smooth(mean(gama{1,1}(1:300,1:199).*active2(1:300,1:199),'omitnan'),20));
hold on;
plot((0.1:0.1:19.9),smooth(mean(gama{2,1}(1:300,1:199).*active2(1:300,1:199),'omitnan'),20));
hold on;
plot((0.1:0.1:19.9),smooth(mean(gama{3,1}(1:300,1:199).*active2(1:300,1:199),'omitnan'),20));
hold on;
plot((0.1:0.1:19.9),smooth(mean(gama{4,1}(1:300,1:199).*active2(1:300,1:199),'omitnan'),20));
hold on;
plot((0.1:0.1:19.9),smooth(mean(gama{5,1}(1:300,1:199).*active2(1:300,1:199),'omitnan'),20));
hold on;
plot([0 20],[0 0],'--','color',[0 0 0]+0.7);
%annotation('textbox',[.60 .68 .1 .1],'String','Flood Tidal Asymmetry','Edgecolor','none','Fontsize',6)
%annotation('textbox',[.60 .32 .1 .1],'String','Ebb Tidal Asymmetry','Edgecolor','none','Fontsize',6)
ylim([-1 1]);set(gca,'Fontname','Times','Fontsize',7);
ylabel('Tidal asymmetry');
xlabel('Distance to upstream boundary (km)')
annotation('textbox',[0.38 0.46 0.06 0.05],'String',{'(d)'},'EdgeColor','none','Fontname','Times');


figure
plot(50*(-50:50),Z1{280}(100:200,180)); hold on;
plot(50*(-50:50),Wl1{280}(100:200,180)-2*P1{2,1}(100:200,180)-2*P1{2,2}(100:200,180))
scatter(50*(-50:50),gama{1,1}(101:201,181))

subchannel1 = gama{1,1}(2:301,2:200).*(gama{1,1}(2:301,2:200)<0);
subchannel1(subchannel1 == 0) = nan ; 
subchannel2 = gama{2,1}(2:301,2:200).*(gama{2,1}(2:301,2:200)<0);
subchannel2(subchannel2 == 0) = nan ; 
subchannel3 = gama{3,1}(2:301,2:200).*(gama{3,1}(2:301,2:200)<0);
subchannel3(subchannel3 == 0) = nan ; 
subchannel4 = gama{4,1}(2:301,2:200).*(gama{4,1}(2:301,2:200)<0);
subchannel4(subchannel4 == 0) = nan ; 
subchannel5 = gama{5,1}(2:301,2:200).*(gama{5,1}(2:301,2:200)<0);
subchannel5(subchannel5 == 0) = nan ; 

figure
plot((0.1:0.1:19.9),smooth(mean(subchannel1,'omitnan'),20));
hold on;
plot((0.1:0.1:19.9),smooth(mean(subchannel2,'omitnan'),20));
hold on;
plot((0.1:0.1:19.9),smooth(mean(subchannel3,'omitnan'),20));
hold on;
plot((0.1:0.1:19.9),smooth(mean(subchannel4,'omitnan'),20));
hold on;
plot((0.1:0.1:19.9),smooth(mean(subchannel5,'omitnan'),20));
hold on;
plot([0 20],[0 0],'--','color',[0 0 0]+0.7);
legend('Base','Scenario 1','Scenario 2','Scenario 3','Scenario 4','location','Southwest');grid on;


print('-dtiff','-r300','Fig7_Tidal_analysis')

%% 
figure('color','white')
set(0,'DefaultAxesFontName','Calibri','DefaultAxesFontSize',10,...
    'DefaultLineLineWidth',1);
set(gcf,'Units','centimeters','Position',[0 0 15.0 10.0]);
colormap(colormap_cpt('balance'));
hold on;
subplot('position',[.1 .7 .3 .15])
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,gama{1,1}(100:200,1:200).*active1(100:200,1:200),100,'linecolor','none');caxis([-1 1]);daspect([1 1 1]);
ylabel('Base')
subplot('position',[.1 .55 .3 .15])
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,gama{2,1}(100:200,1:200).*active1(100:200,1:200),100,'linecolor','none');caxis([-1 1]);daspect([1 1 1]);
ylabel('Sce 1')
subplot('position',[.1 .4 .3 .15])
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,gama{3,1}(100:200,1:200).*active1(100:200,1:200),100,'linecolor','none');caxis([-1 1]);daspect([1 1 1]);
ylabel('Sce 2')
subplot('position',[.1 .25 .3 .15])
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,gama{4,1}(100:200,1:200).*active1(100:200,1:200),100,'linecolor','none');caxis([-1 1]);daspect([1 1 1]);
ylabel('Sce 3')
subplot('position',[.1 .1 .3 .15])
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,gama{5,1}(100:200,1:200).*active1(100:200,1:200),100,'linecolor','none');caxis([-1 1]);daspect([1 1 1]);
ylabel('Sce 4')
xlabel('Distance (km)')
colorbar('position',[.43 .12 .02 .72])
subplot('position',[.55 .3 .4 .5])
plot((0.1:0.1:20),smooth(mean(gama{1,1}(1:300,1:200).*active2(1:300,1:200),'omitnan'),20));
hold on;
plot((0.1:0.1:20),smooth(mean(gama{2,1}(1:300,1:200).*active2(1:300,1:200),'omitnan'),20));
hold on;
plot((0.1:0.1:20),smooth(mean(gama{3,1}(1:300,1:200).*active2(1:300,1:200),'omitnan'),20));
hold on;
plot((0.1:0.1:20),smooth(mean(gama{4,1}(1:300,1:200).*active2(1:300,1:200),'omitnan'),20));
hold on;
plot((0.1:0.1:20),smooth(mean(gama{5,1}(1:300,1:200).*active2(1:300,1:200),'omitnan'),20));
hold on;
plot([0 20],[0 0],'--','color',[0 0 0]+0.2);
annotation('textbox',[.55 .68 .1 .1],'String','Flood Tidal Asymmetry','Edgecolor','none','Fontsize',8)
annotation('textbox',[.55 .32 .1 .1],'String','Ebb Tidal Asymmetry','Edgecolor','none','Fontsize',8)
ylim([-1 1]);legend({'Base model','Scenario 1','Scenario 2','Scenario 3','Scenario 4'},'location','SouthEast');
title('Tidal asymmetry');
xlabel('Distance to upstream boundary (km)')
grid on
print fig8.jpg -djpeg -r600


%% Fig 9���ӣţģɣͣţΣԡ��Уңϣأ�
figure('color','white')
set(gcf,'Units','centimeters','Position',[0 0 19.0 15.0]);
subplot(2,1,1)
plot((0.1:0.1:20),smooth(sum(P1{4,1}(1:300,1:200),'omitnan'),50),'linewidth',1.5,'color','#0072BD');hold on;
plot((0.1:0.1:20),smooth(sum(P2{4,1}(1:300,1:200),'omitnan'),50),'linewidth',1.5,'color','#D95319');
plot((0.1:0.1:20),smooth(sum(P3{4,1}(1:300,1:200),'omitnan'),50),'linewidth',1.5,'color','#EDB120');
plot((0.1:0.1:20),smooth(sum(P4{4,1}(1:300,1:200),'omitnan'),50),'linewidth',1.5,'color','#7E2F8E');
plot((0.1:0.1:20),smooth(sum(P5{4,1}(1:300,1:200),'omitnan'),50),'linewidth',1.5,'color','#77AC30');
plot((0.1:0.1:20),smooth(sum(P1{5,1}(1:300,1:200),'omitnan'),50),'-.','linewidth',1.5,'color','#0072BD');
plot((0.1:0.1:20),smooth(sum(P2{5,1}(1:300,1:200),'omitnan'),50),'-.','linewidth',1.5,'color','#D95319');
plot((0.1:0.1:20),smooth(sum(P3{5,1}(1:300,1:200),'omitnan'),50),'-.','linewidth',1.5,'color','#EDB120');
plot((0.1:0.1:20),smooth(sum(P4{5,1}(1:300,1:200),'omitnan'),50),'-.','linewidth',1.5,'color','#7E2F8E');
plot((0.1:0.1:20),smooth(sum(P5{5,1}(1:300,1:200),'omitnan'),50),'-.','linewidth',1.5,'color','#77AC30');
plot((0.1:0.1:20),smooth(sum(P1{6,1}(1:300,1:200),'omitnan'),50),'--','linewidth',1.5,'color','#0072BD');
plot((0.1:0.1:20),smooth(sum(P2{6,1}(1:300,1:200),'omitnan'),50),'--','linewidth',1.5,'color','#D95319');
plot((0.1:0.1:20),smooth(sum(P3{6,1}(1:300,1:200),'omitnan'),50),'--','linewidth',1.5,'color','#EDB120');
plot((0.1:0.1:20),smooth(sum(P4{6,1}(1:300,1:200),'omitnan'),50),'--','linewidth',1.5,'color','#7E2F8E');
plot((0.1:0.1:20),smooth(sum(P5{6,1}(1:300,1:200),'omitnan'),50),'--','linewidth',1.5,'color','#77AC30');grid on;ylim([-2 5]);
%'Q_{ent},Eulerian residual current','Q_{eia},Euler-induced asymmetry','Q_{tia},Tidal-induced asymmetry',
%scatter(index_tq1/10,4,[],'v','filled','MarkerEdgeColor','#0072BD','MarkerFaceColor','#0072BD');hold on;
%scatter(index_tq2/10,4,[],'v','filled','MarkerEdgeColor','#D95319','MarkerFaceColor','#D95319');
%scatter(index_tq3/10,4,[],'v','filled','MarkerEdgeColor','#EDB120','MarkerFaceColor','#EDB120');
%scatter(index_tq4/10,4,[],'v','filled','MarkerEdgeColor','#7E2F8E','MarkerFaceColor','#7E2F8E');
%scatter(index_tq5/10,4,[],'v','filled','MarkerEdgeColor','#77AC30','MarkerFaceColor','#77AC30');
set(gca,'Fontname','Times','Fontsize',10); 
annotation('textbox',[0.15 0.87 0.06 0.05],'String',{'(a)'},'EdgeColor','none','Fontname','Times');
subplot(2,1,2)
%plot([index_tq1/10 index_tq1e/10],[4.5 4.5],'color','#0072BD');
%plot([index_tq2/10 index_tq2e/10],[4.75 4.75],'color','#D95319');
%plot([index_tq3/10 index_tq3e/10],[5 5],'color','#EDB120');
%plot([index_tq4/10 index_tq4e/10],[5.25 5.25],'color','#7E2F8E');
%plot([index_tq5/10 index_tq5e/10],[5.5 5.5],'color','#77AC30');
%total residual sediment transport
h1 = plot((0.1:0.1:20),smooth(sum((sum_sed_prx1),'omitnan'),50),'linewidth',1.5,'color','#0072BD');hold on;
h2 = plot((0.1:0.1:20),smooth(sum((sum_sed_prx2),'omitnan'),50),'linewidth',1.5,'color','#D95319');
h3 = plot((0.1:0.1:20),smooth(sum((sum_sed_prx3),'omitnan'),50),'linewidth',1.5,'color','#EDB120');
h4 = plot((0.1:0.1:20),smooth(sum((sum_sed_prx4),'omitnan'),50),'linewidth',1.5,'color','#7E2F8E');
h5 = plot((0.1:0.1:20),smooth(sum((sum_sed_prx5),'omitnan'),50),'linewidth',1.5,'color','#77AC30');
%ylim([-2 5]);
set(gca,'Fontname','Times','Fontsize',10); xlabel('Distance to upstream boundary (km)')
yb = ylabel('Dimensionless Sediement Transport Rate (-)');
yb.Position(2) = yb.Position(2)+1.7;
grid on; annotation('textbox',[0.15 0.40 0.06 0.05],'String',{'(b)'},'EdgeColor','none','Fontname','Times');

legend([h1 h2 h3 h4 h5],{'Base model','Scenario 1','Scenario 2','Scenario 3','Scenario 4'},'location','Southeast');
print fig9.jpg -djpeg -r600

figure('color','white')
set(0,'DefaultAxesFontName','Calibri','DefaultAxesFontSize',7,...
    'DefaultLineLineWidth',1);
set(gcf,'Units','centimeters','Position',[0 0 15.0 9.0]);
colormap(colormap_cpt('seismic'));
k = [-0.5 0.5];
subplot(5,3,1)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P1{4,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);
ylabel('Base Model'); set(gca,'xtick',[])
subplot(5,3,4)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P2{4,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);
ylabel('Scenario 1');set(gca,'xtick',[])
subplot(5,3,7)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P3{4,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);
ylabel('Scenario 2');set(gca,'xtick',[])
subplot(5,3,10)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P4{4,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);
ylabel('Scenario 3');set(gca,'xtick',[])
subplot(5,3,13)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P5{4,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);
ylabel('Scenario 4');
xlabel('Q_{ent}');
subplot(5,3,2)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P1{5,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'xtick',[])
subplot(5,3,5)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P2{5,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'xtick',[])
subplot(5,3,8)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P3{5,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'xtick',[])
subplot(5,3,11)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P4{5,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'xtick',[])
subplot(5,3,14)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P5{5,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);
xlabel('Q_{eia}');
hold on;
subplot(5,3,3)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P1{6,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'xtick',[])
subplot(5,3,6)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P2{6,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'xtick',[])
subplot(5,3,9)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P3{6,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'xtick',[])
subplot(5,3,12)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P4{6,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);set(gca,'xtick',[])
subplot(5,3,15)
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,P5{6,1}(100:200,1:200),100,'linecolor','none');caxis(k);daspect([1 1 1]);
colorbar('position',[.92 .12 .02 .8]);
xlabel('Q_{tia}');
print fig91.jpg -djpeg -r600


figure('color','white')
set(0,'DefaultAxesFontName','Calibri','DefaultAxesFontSize',7,...
    'DefaultLineLineWidth',1);
set(gcf,'Units','centimeters','Position',[0 0 15.0 10.0]);
colormap(colormap_cpt('balance'));
hold on;
subplot('position',[.1 .7 .3 .15])
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,gama{1,1}(100:200,1:200).*active1(100:200,1:200),100,'linecolor','none');caxis([-1 1]);daspect([1 1 1]);
subplot('position',[.1 .55 .3 .15])
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,gama{2,1}(100:200,1:200).*active1(100:200,1:200),100,'linecolor','none');caxis([-1 1]);daspect([1 1 1]);
subplot('position',[.1 .4 .3 .15])
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,gama{3,1}(100:200,1:200).*active1(100:200,1:200),100,'linecolor','none');caxis([-1 1]);daspect([1 1 1]);
subplot('position',[.1 .25 .3 .15])
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,gama{4,1}(100:200,1:200).*active1(100:200,1:200),100,'linecolor','none');caxis([-1 1]);daspect([1 1 1]);
subplot('position',[.1 .1 .3 .15])
contourf(X(100:200,1:200)/1000,Y(100:200,1:200)/1000,gama{5,1}(100:200,1:200).*active1(100:200,1:200),100,'linecolor','none');caxis([-1 1]);daspect([1 1 1]);
colorbar('position',[.43 .12 .02 .72])
subplot('position',[.55 .3 .4 .5])
plot((0.1:0.1:20),smooth(mean(gama{1,1}(1:300,1:200).*active2(1:300,1:200),'omitnan'),20));
hold on;
plot((0.1:0.1:20),smooth(mean(gama{2,1}(1:300,1:200).*active2(1:300,1:200),'omitnan'),20));
hold on;
plot((0.1:0.1:20),smooth(mean(gama{3,1}(1:300,1:200).*active2(1:300,1:200),'omitnan'),20));
hold on;
plot((0.1:0.1:20),smooth(mean(gama{4,1}(1:300,1:200).*active2(1:300,1:200),'omitnan'),20));
hold on;
plot((0.1:0.1:20),smooth(mean(gama{5,1}(1:300,1:200).*active2(1:300,1:200),'omitnan'),20));
hold on;
plot([0 20],[0 0],'--','color',[0 0 0]+0.2);
annotation('textbox',[.60 .68 .1 .1],'String','Flood Tidal Asymmetry','Edgecolor','none','Fontsize',6)
annotation('textbox',[.60 .32 .1 .1],'String','Ebb Tidal Asymmetry','Edgecolor','none','Fontsize',6)
ylim([-1 1]);legend({'Base model','Scenario 1','Scenario 2','Scenario 3','Scenario 4'},'location','SouthEast');
title('Tidal asymmetry');
xlabel('Distance to upstream boundary (km)')
grid on
print('-dtiff','-r600',[figsavedir 'tidal asymmetry (gama)'])


%% Fig 9 cumulative sediment transport
figure('color','white')
set(gcf,'Units','centimeters','Position',[0 0 7 5.0]);
plot(0:100/10512:100,movmean(cell2mat(SBTRC_sand2{6})-cell2mat(SBTRC_sand1{6}),149),'linewidth',1.5);hold on;
plot(0:100/10512:100,movmean(cell2mat(SBTRC_sand3{6})-cell2mat(SBTRC_sand1{6}),149),'linewidth',1.5);
plot(0:100/10512:100,movmean(cell2mat(SBTRC_sand4{6})-cell2mat(SBTRC_sand1{6}),149),'linewidth',1.5);
plot(0:100/10512:100,movmean(cell2mat(SBTRC_sand5{6})-cell2mat(SBTRC_sand1{6}),149),'linewidth',1.5);
xlabel('Time (year)');ylabel('Cumulative sediment transport (m^3)');set(gca,'Fontname','Times','Fontsize',7);
legend('Scenario 1','Scenario 2','Scenario 3','Scenario 4','location','Southwest','box','off')
print fig10.jpg -djpeg -r600


%% Fig 10 velocity_decompose
figure('color','white')
set(0,'DefaultAxesFontName','Calibri','DefaultAxesFontSize',7,...
    'DefaultLineLineWidth',1);
set(gcf,'Units','centimeters','Position',[0 0 15.0 18.0]);

x= 0.13; y = 0.905; dx = 0.2; dy = 0.22;
annotation('textbox',[x y 0.06 0.05],'String',{'(a)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x y-dy 0.06 0.05],'String',{'(b)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x y-dy*2 0.06 0.05],'String',{'(c)'},'EdgeColor','none','Fontname','Times','Fontsize',8);
annotation('textbox',[x y-dy*3 0.06 0.05],'String',{'(d)'},'EdgeColor','none','Fontname','Times','Fontsize',8);

subplot(4,1,1);hold on; k = [-0.04 0.04];
title('scenario 1')
bar(5:5:95,[VD_c2.*sqrt(VD_r1.*VD_s1)-VD_c1.*sqrt(VD_r1.*VD_s1) 
          VD_c1.*sqrt(VD_r2.*VD_s1)-VD_c1.*sqrt(VD_r1.*VD_s1) 
          VD_c1.*sqrt(VD_r1.*VD_s2)-VD_c1.*sqrt(VD_r1.*VD_s1)]);
plot((5:5:95),VD_c2.*sqrt(VD_r2.*VD_s2)-VD_c1.*sqrt(VD_r1.*VD_s1),'-o','markersize',4,'linewidth',1.5);grid on;ylim(k)
plot((5:5:95),VD_v2-VD_v1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'Fontname','Times','Fontsize',7);ylabel('velocity difference (m/s)')
subplot(4,1,2);hold on; ylabel('velocity difference (m/s)')
title('scenario 2')
bar(5:5:95,[VD_c3.*sqrt(VD_r1.*VD_s1)-VD_c1.*sqrt(VD_r1.*VD_s1) 
          VD_c1.*sqrt(VD_r3.*VD_s1)-VD_c1.*sqrt(VD_r1.*VD_s1) 
          VD_c1.*sqrt(VD_r1.*VD_s3)-VD_c1.*sqrt(VD_r1.*VD_s1)]);
plot((5:5:95),VD_c3.*sqrt(VD_r3.*VD_s3)-VD_c1.*sqrt(VD_r1.*VD_s1),'-o','markersize',4,'linewidth',1.5);grid on;ylim(k)
plot((5:5:95),VD_v3-VD_v1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'Fontname','Times','Fontsize',7);
subplot(4,1,3);hold on;ylabel('velocity difference (m/s)')
title('scenario 3')
bar(5:5:95,[VD_c4.*sqrt(VD_r1.*VD_s1)-VD_c1.*sqrt(VD_r1.*VD_s1) 
          VD_c1.*sqrt(VD_r4.*VD_s1)-VD_c1.*sqrt(VD_r1.*VD_s1) 
          VD_c1.*sqrt(VD_r1.*VD_s4)-VD_c1.*sqrt(VD_r1.*VD_s1)]);
plot((5:5:95),VD_c4.*sqrt(VD_r4.*VD_s4)-VD_c1.*sqrt(VD_r1.*VD_s1),'-o','markersize',4,'linewidth',1.5);grid on;ylim(k)
plot((5:5:95),VD_v4-VD_v1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'Fontname','Times','Fontsize',7);
subplot(4,1,4);hold on;ylabel('velocity difference (m/s)')
title('scenario 4')
bar(5:5:95,[VD_c5.*sqrt(VD_r1.*VD_s1)-VD_c1.*sqrt(VD_r1.*VD_s1) 
          VD_c1.*sqrt(VD_r5.*VD_s1)-VD_c1.*sqrt(VD_r1.*VD_s1) 
          VD_c1.*sqrt(VD_r1.*VD_s5)-VD_c1.*sqrt(VD_r1.*VD_s1)]);
plot((5:5:95),VD_c5.*sqrt(VD_r5.*VD_s5)-VD_c1.*sqrt(VD_r1.*VD_s1),'-o','markersize',4,'linewidth',1.5);grid on;ylim(k)
plot((5:5:95),VD_v5-VD_v1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'Fontname','Times','Fontsize',7);
legend('Chezy Roughness','Water Depth','Surface Slope','Total difference','Model result','position',[0.26,0.87,0.55,0.05],'NumColumns',3,'Box','off')
xlabel('Time (year)');ylabel('velocity difference (m/s)')
print fig11.eps -depsc2 -r600


figure('color','white')
set(0,'DefaultAxesFontName','Calibri','DefaultAxesFontSize',7,...
    'DefaultLineLineWidth',1);
set(gcf,'Units','centimeters','Position',[0 0 15.0 18.0]);
subplot(4,1,1);hold on; k = [-0.04 0.04];
title('scenario 1')
bar(5:5:95,[VD_m2.^-1.*VD_r1.^(2/3).*VD_s1.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5 
          VD_m1.^-1.*VD_r2.^(2/3).*VD_s1.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5 
          VD_m1.^-1.*VD_r1.^(2/3).*VD_s2.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5]);
plot((5:5:95),VD_m2.^-1.*VD_r2.^(2/3).*VD_s2.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5,'-o','markersize',4,'linewidth',1.5);grid on;ylim(k)
plot((5:5:95),VD_v2-VD_v1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'Fontname','Times','Fontsize',7);ylabel('velocity difference (m/s)')
subplot(4,1,2);hold on; ylabel('velocity difference (m/s)')
title('scenario 2')
bar(5:5:95,[VD_m3.^-1.*VD_r1.^(2/3).*VD_s1.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5 
          VD_m1.^-1.*VD_r3.^(2/3).*VD_s1.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5 
          VD_m1.^-1.*VD_r1.^(2/3).*VD_s3.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5]);
plot((5:5:95),VD_m3.^-1.*VD_r3.^(2/3).*VD_s3.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5,'-o','markersize',4,'linewidth',1.5);grid on;ylim(k)
plot((5:5:95),VD_v3-VD_v1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'Fontname','Times','Fontsize',7);
subplot(4,1,3);hold on;ylabel('velocity difference (m/s)')
title('scenario 3')
bar(5:5:95,[VD_m4.^-1.*VD_r1.^(2/3).*VD_s1.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5 
          VD_m1.^-1.*VD_r4.^(2/3).*VD_s1.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5 
          VD_m1.^-1.*VD_r1.^(2/3).*VD_s4.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5]);
plot((5:5:95),VD_m4.^-1.*VD_r4.^(2/3).*VD_s4.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5,'-o','markersize',4,'linewidth',1.5);grid on;ylim(k)
plot((5:5:95),VD_v4-VD_v1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'Fontname','Times','Fontsize',7);
subplot(4,1,4);hold on;ylabel('velocity difference (m/s)')
title('scenario 4')
bar(5:5:95,[VD_m5.^-1.*VD_r1.^(2/3).*VD_s1.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5 
          VD_m1.^-1.*VD_r5.^(2/3).*VD_s1.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5 
          VD_m1.^-1.*VD_r1.^(2/3).*VD_s5.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5]);
plot((5:5:95),VD_m5.^-1.*VD_r5.^(2/3).*VD_s5.^0.5-VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5,'-o','markersize',4,'linewidth',1.5);grid on;ylim(k)
plot((5:5:95),VD_v5-VD_v1,'--','linewidth',1.5,'color',[0.2 0.2 0.2]);set(gca,'Fontname','Times','Fontsize',7);
legend("Mannings Roughness",'Water Depth','Surface Slope','Total difference','Model result','position',[0.26,0.87,0.55,0.05],'NumColumns',3,'Box','off')
xlabel('Time (year)');ylabel('velocity difference (m/s)')
print('-dpng','-r300','fig10_velocity_decompose')


figure
plot(VD_m5.^-1.*VD_r5.^(2/3).*VD_s5.^0.5); hold on
plot(VD_v5)

figure
plot(VD_m1.^-1.*VD_r1.^(2/3).*VD_s1.^0.5); hold on
plot(VD_v1)


%% flow resistance

figure('color','white')
set(0,'DefaultAxesFontName','Calibri','DefaultAxesFontSize',14,...
    'DefaultLineLineWidth',1);
set(gcf,'Units','centimeters','Position',[0 0 15.0 6.0]); 
subplot(1,5,5)
plot((5:5:95),VD_r5.^(1/6)./VD_c5,'-o');hold on;
plot((5:5:95),VD_r5.^(1/6)./(3.1*VD_r5.^0.33./(VD_s5.^0.38)),'-o');
plot((5:5:95),VD_r5.^(1/6)./(7.14*VD_r5.^0.17./(VD_s5.^0.17)),'-o');
plot((5:5:95),VD_r5.^(1/6)./(6.04*VD_r5.^0.32./(VD_s5.^0.24)),'-o');ylim([0 0.04])
legend('Model result','Jarrett (1984)','Bjerklie et al.(2005)','Lopez et al.(2007)'); 
subplot(1,5,1)
plot((5:5:95),VD_r1.^(1/6)./VD_c1,'-o');hold on;
plot((5:5:95),VD_r1.^(1/6)./(3.1*VD_r1.^0.33./(VD_s1.^0.38)),'-o');
plot((5:5:95),VD_r1.^(1/6)./(7.14*VD_r1.^0.17./(VD_s1.^0.17)),'-o');
plot((5:5:95),VD_r1.^(1/6)./(6.04*VD_r1.^0.32./(VD_s1.^0.24)),'-o');ylim([0 0.04])
subplot(1,5,2)
plot((5:5:95),VD_r2.^(1/6)./VD_c2,'-o');hold on;
plot((5:5:95),VD_r2.^(1/6)./(3.1*VD_r2.^0.33./(VD_s2.^0.38)),'-o');
plot((5:5:95),VD_r2.^(1/6)./(7.14*VD_r2.^0.17./(VD_s2.^0.17)),'-o');
plot((5:5:95),VD_r2.^(1/6)./(6.04*VD_r2.^0.32./(VD_s2.^0.24)),'-o');ylim([0 0.04])
subplot(1,5,3)
plot((5:5:95),VD_r3.^(1/6)./VD_c3,'-o');hold on;
plot((5:5:95),VD_r3.^(1/6)./(3.1*VD_r3.^0.33./(VD_s3.^0.38)),'-o');
plot((5:5:95),VD_r3.^(1/6)./(7.14*VD_r3.^0.17./(VD_s3.^0.17)),'-o');
plot((5:5:95),VD_r3.^(1/6)./(6.04*VD_r3.^0.32./(VD_s3.^0.24)),'-o');ylim([0 0.04])
subplot(1,5,4)
plot((5:5:95),VD_r4.^(1/6)./VD_c4,'-o');hold on;
plot((5:5:95),VD_r4.^(1/6)./(3.1*VD_r4.^0.33./(VD_s4.^0.38)),'-o');
plot((5:5:95),VD_r4.^(1/6)./(7.14*VD_r4.^0.17./(VD_s4.^0.17)),'-o');
plot((5:5:95),VD_r4.^(1/6)./(6.04*VD_r4.^0.32./(VD_s4.^0.24)),'-o');ylim([0 0.04])

v = [VD_v1 VD_v2 VD_v3 VD_v4 VD_v5]
r = [VD_r1 VD_r2 VD_r3 VD_r4 VD_r5 ]
s = [VD_s1 VD_s2 VD_s3 VD_s4 VD_s5]



% Morphological response time
for i = 1:70
    Prism01(i) = sum(Q1{19}(i*149-148:i*149).*(Q1{19}(i*149-148:i*149)<0))*-5*60;
    Prism02(i) = sum(Q2{19}(i*149-148:i*149).*(Q2{19}(i*149-148:i*149)<0))*-5*60;
    Prism03(i) = sum(Q3{19}(i*149-148:i*149).*(Q3{19}(i*149-148:i*149)<0))*-5*60;
    Prism04(i) = sum(Q4{19}(i*149-148:i*149).*(Q4{19}(i*149-148:i*149)<0))*-5*60;
    Prism05(i) = sum(Q5{19}(i*149-148:i*149).*(Q5{19}(i*149-148:i*149)<0))*-5*60;

end
figure
plot(Prism01); hold on;
plot(Prism02);
plot(Prism03);
plot(Prism04);
plot(Prism05);
xlabel('Time (number of tidal cycles)');
ylabel('Tidal Prism (m^3)')
legend('Base','Scenario 1','Scenario 2','Scenario 3','Scenario 4','location','Southwest');grid on;


    area01 = sum(sum(Z1{283}(1:300,1:199)<0))*5000;
    ve01 = 65*10^-6*(Prism01./10^6).^(3/2)*10^6;
    tao01 = 1/2/(7.5*10^-6)*(ve01./10);
    
    figure
    plot(tao01/3600/24/365)
    xlabel('Time (number of tidal cycles)');
    ylabel('Tidal response time(yr)')

